Starting point¶
Maximum tidal amplitude compared against TPXO 7.2 solution¶
What I try do achive is a comparison of maximum tidal amplitude from ROMS against the TPXO7.2 solution (provides forcing for ROMS). To do so, I compare ROMS maximum sea surface height anomaly from a hourly record of 14 days against the summation of the amplitudes over all (8) constituents from the forcing file.
import xarray as xr
import matplotlib.pyplot as plt
import os
import sys
import numpy as np
from dotenv import load_dotenv, find_dotenv
# find .env automagically by walking up directories until it's found
dotenv_path = find_dotenv()
# load up the entries as environment variables
load_dotenv(dotenv_path)
src_dir = os.environ.get('srcdir')
print(src_dir)
sys.path.append(src_dir)
# always reload modules marked with "%aimport"
%load_ext autoreload
%autoreload 1
%aimport features.grid_ttide
from features.grid_ttide import grid_ttide, plot_M2O1_diff
file_path = os.path.join(os.environ.get('rawdir'),'waom10_full_forcing','ocean_his_hourly.nc')
ds_ref = xr.open_mfdataset(file_path,chunks={'eta_rho':100,'xi_rho':100})
# calculate maximum anomaly of the free surface height compared to the 14 day mean
ref_maxAmpl = ds_ref.zeta.max('ocean_time')-ds_ref.zeta.mean('ocean_time')
# load tpxo forcing amplitude and phase information as Xarray dataset
tpxo_path = os.path.join(os.environ.get('rawdir'),'waom10_full_forcing','Data','waom10','waom10_tds_large.nc')
dsf = xr.open_mfdataset(tpxo_path)
# calculate the superposition of amplitudes of all constituents (max tidal amplitude possible)
tpxo_maxAmpl = dsf.sum('tide_period')
# plot both next to each other with same colorbounds
plt.clf()
fig,(ax1,ax2) = plt.subplots(ncols=2,figsize=(15,5))
fig.suptitle('Comparison of maximum tidal height amplitude ROMS vs. TPXO7.2',fontsize=16)
ref_maxAmpl.fillna(0).plot(ax=ax1,vmax=4)
ax1.set_title('ROMS (from 14 day record)')
ax1.axis('off')
tpxo_maxAmpl.tide_Eamp.plot(ax=ax2,vmax=4)
ax2.set_title('TPXO (theoretical)')
ax2.axis('off')
plt.show()
... Roms gets about two times to high under the big ice shelves. Since the height looks ok everywhere else, this might be related to the kelvin waves. Two ideas to fix this:
- no forcing on the continental shelf
- increase bottom/ice shelf drag, high resolution might help
Comparison of M2 and O1 amplitude maps¶
(others violate Rayleigh condition with just 14 day record) Apply t_tide tidal harmonic analysis on every grid point and fill O1 and M2 amplitude maps in Xarray dataset. First define the function
ta.grid_ttide(ds_ref,ds_ref)
then apply to dataset.
grid_lr_path = os.path.join(os.environ.get('rawdir'),'waom10_full_forcing','Data','waom10','waom10_grd_large.nc')
grid_lr = xr.open_mfdataset(grid_lr_path)
ds_ref = ta.grid_ttide(ds_ref,grid_lr,res=50)
Plot against TPXO O1 and M2 amplitude maps
plt.close('all')
fig,axes = plt.subplots(2,2,figsize=(15,12))
ax1,ax2,ax3,ax4 = axes.flatten()
fig.suptitle('Comparison of M2 and O1 height amplitude ROMS vs. TPXO7.2',fontsize=16)
ds_ref.M2_amp.fillna(0).plot(ax=ax1,vmax=1.5)
ax1.set_title('ROMS M2 amplitude in m')
ax1.axis('off')
dsf.tide_Eamp[0].plot(ax=ax2,vmax=1.5)
ax2.set_title('TPXO M2 amplitude in m')
ax2.axis('off')
ds_ref.O1_amp.fillna(0).plot(ax=ax3,vmax=0.5)
ax3.set_title('ROMS O1 amplitude in m')
ax3.axis('off')
dsf.tide_Eamp[5].plot(ax=ax4,vmax=0.5)
ax4.set_title('TPXO O1 amplitude in m')
ax4.axis('off')
plt.show()
... both are unacceptable to high.
#mask tidal forcing
pott_path = os.path.join(os.environ.get('rawdir'),'waom10_full_forcing','Data','waom10','waom10_ptds_large.nc')
pott_ds = xr.open_dataset(pott_path)
grid_path = os.path.join(os.environ.get('rawdir'),'waom10_full_forcing','Data','waom10','waom10_grd_large.nc')
grid_ds = xr.open_dataset(grid_path)
Pamp_masked = pott_ds.tide_Pamp.where(grid_lr.zice.values == 0.0 ,0.0)
plt.close()
Pamp_masked[2].plot()
plt.show()
Effect on max surface height anomaly¶
plot against tpxo as before, but define function for later plots
#look at results
file_path = os.path.join(os.environ.get('rawdir'),'waom10_tidal_mask','ocean_his_tm_0006.nc')
ds_tm = xr.open_mfdataset(file_path,chunks={'eta_rho':100,'xi_rho':100})
# calculate maximum anomaly of the free surface height compared to the 14 day mean
tm_maxAmpl = ds_tm.zeta[336:].max('ocean_time')-ds_tm.zeta[336:].mean('ocean_time')
def plot_max_ampl(case_max_ampl_da,case_str,ref_max_ampl,ref_str):
plt.clf()
fig,(ax1,ax2,ax3) = plt.subplots(ncols=3,figsize=(17,4))
fig.suptitle('Comparison of maximum tidal height amplitude ROMS vs. TPXO7.2',fontsize=16)
ref_max_ampl.fillna(0).plot(ax=ax1,vmax=4)
ax1.set_title(ref_str+' [m]')
ax1.axis('off')
case_max_ampl_da.fillna(0).plot(ax=ax2,vmax=4)
ax2.set_title(case_str+ ' [m]')
ax2.axis('off')
tpxo_maxAmpl.tide_Eamp.plot(ax=ax3,vmax=4)
ax3.set_title('TPXO (theoretical)')
ax3.axis('off')
plt.show()
plot_max_ampl(tm_maxAmpl,'tidal mask',ref_maxAmpl,'base case')
... this looks a lot better, though its still ssh amplitudes are still to high in the shallow areas under the ice shelves.
M2 and O1 maps¶
ds_tm = grid_ttide(ds_tm,grid_lr,50)
# load tpxo forcing amplitude and phase information as Xarray dataset
def plot_M2O1(case_ds,case_str,ds_ref,ref_str):
plt.close('all')
fig,axes = plt.subplots(ncols=3,nrows=2,figsize=(17,8))
ax1,ax2,ax3,ax4,ax5,ax6 = axes.flatten()
fig.suptitle('Comparison of M2 and O1 height amplitude',fontsize=16)
ds_ref.M2_ampl.fillna(0).plot(ax=ax1,vmax=1.5)
ax1.set_title(ref_str+' [m]')
ax1.axis('off')
case_ds.M2_ampl.fillna(0).plot(ax=ax2,vmax=1.5)
ax2.set_title(case_str +' [m]')
ax2.axis('off')
dsf.tide_Eamp[0].plot(ax=ax3,vmax=1.5)
ax3.set_title('TPXO M2 amplitude in m')
ax3.axis('off')
ds_ref.O1_ampl.fillna(0).plot(ax=ax4,vmax=0.5)
ax4.set_title(ref_str+' [m]')
ax4.axis('off')
case_ds.O1_ampl.fillna(0).plot(ax=ax5,vmax=0.5)
ax5.set_title(case_str+' [m]')
ax5.axis('off')
dsf.tide_Eamp[5].plot(ax=ax6,vmax=0.5)
ax6.set_title('TPXO O1 amplitude in m')
ax6.axis('off')
plt.show()
plot_M2O1(ds_tm,'tidal mask',ds_ref,'base case')
... seems like we get a big step at the mask boundary.
5x Drag coefficient¶
Ramping up the drag coefficent from 0.003 to 0.015 to represent conversion to internal tides, which is not present in 10 km model run.
#look at results
file_path = os.path.join(os.environ.get('rawdir'),'waom10_tidal_mask','ocean_his_5xCd_0006.nc')
ds_Cd = xr.open_dataset(file_path,chunks={'eta_rho':100,'xi_rho':100})
# calculate maximum anomaly of the free surface height compared to the 14 day mean
Cd_maxAmpl = ds_Cd.zeta[336:].max('ocean_time')-ds_Cd.zeta[336:].mean('ocean_time')
plot_max_ampl(Cd_maxAmpl,'5xCd',ref_maxAmpl,'base case')
... visible, but small decrease in max amplitude.
2x viscosity¶
Ramping up the visc2 from 500 m^2/s to 1000 m^2/s.
#look at results
import os
import xarray as xr
file_path = os.path.join(os.environ.get('rawdir'),'waom10_tidal_mask','ocean_his_2xVisc2_0006.nc')
ds_Vi = xr.open_dataset(file_path,chunks={'eta_rho':100,'xi_rho':100})
# calculate maximum anomaly of the free surface height compared to the 14 day mean
Vi_maxAmpl = ds_Vi.zeta[336:].max('ocean_time')-ds_Vi.zeta[336:].mean('ocean_time')
plot_max_ampl(Vi_maxAmpl,'2x viscosity',ref_maxAmpl,'base case')
... almost no impact.
Mask shallow areas¶
Since bathymetry has such a big impact on tides, maybe our artificial deepening at the grounding lines create the bias.
Masking shallow areas (including all artificially deepened cells) by setting mask_rho == 0 where wct <= 50 m.
#look at results
import os
import xarray as xr
file_path = os.path.join(os.environ.get('rawdir'),'waom10_tidal_mask','ocean_his_ms_0006.nc')
ds_ms = xr.open_dataset(file_path,chunks={'eta_rho':100,'xi_rho':100})
# calculate maximum anomaly of the free surface height compared to the 14 day mean
ms_maxAmpl = ds_ms.zeta[336:].max('ocean_time')-ds_ms.zeta[336:].mean('ocean_time')
plot_max_ampl(ms_maxAmpl,'mask shallow',ref_maxAmpl,'base case')
... a bit better in the Weddell Sea and Amery. However in the Ross sea it seems just to shift the high amplitudes further in the direction of the kelvin wave propagation. ... more ideas would be trying a different vertical mixing scheme as recommended by Robin Robertson (MY and some GLS) and looking at the impact of increased resolution.
2x horizontal resolution¶
By increasing the horizontal resolution from 10km to 5km we hope to lose energy into baroclinic tides, which we don't resolve withthe coarse grid. Also better representation of bathymetry and ice draft, less smoothing and artificial deepening to only 30m water column thickness might have an impact.
#look at results
file_path = os.path.join(os.environ.get('rawdir'),'gdata','waom5','ocean_his_0002.nc')
ds_hr = xr.open_mfdataset(file_path,chunks={'xi_rho':200,'eta_rho':200})
# calculate maximum anomaly of the free surface height compared to the 14 day mean
hr_maxAmpl = ds_hr.zeta.max('ocean_time')-ds_hr.zeta.mean('ocean_time')
Max tidal height amplitude¶
plot_max_ampl(hr_maxAmpl,'2x horizontal resolution',ref_maxAmpl,'base case')
M2 O1 height amplitudes¶
ds_hr = grid_ttide(ds_hr,ds_hr,50)
plot_M2O1(ds_hr,'2x horizontal resolution',ds_ref,'base case')
2x horizontal resolution and MY25 vertical mixing scheme¶
Following advice from Robin Roberston, I've changed the vertical mixing from LMD to MY. However, we expect MY to mix less than LMD and therefore possibly dissipate less energy?
#look at results
file_path = os.path.join(os.environ.get('rawdir'),'gdata','waom5_MY','ocean_his_0002.nc')
ds_my = xr.open_mfdataset(file_path,chunks={'xi_rho':200,'eta_rho':200})
# calculate maximum anomaly of the free surface height compared to the 14 day mean
my_maxAmpl = ds_my.zeta[:].max('ocean_time')-ds_my.zeta[:].mean('ocean_time')
Max tidal height amplitude¶
plot_max_ampl(my_maxAmpl,'high res + MY25 mixing',hr_maxAmpl,'high res')
M2 O1 height amplitudes¶
ds_my = grid_ttide(ds_my,ds_my,50)
plot_M2O1(ds_my,'2x hr and MY25 mixing',ds_hr,'high res')
2x horizontal resolution and no topography smoothing¶
Following Robins papaers, I try a case, where I don't apply any smooting to the topographz (bed rock and ice draft) to enhance internal tide generation over rough topography.
#look at results
file_path = os.path.join(os.environ.get('rawdir'),'waom5_rough','ocean_his_0002.nc')
ds_ns = xr.open_mfdataset(file_path,chunks={'xi_rho':200,'eta_rho':200})
# calculate maximum anomaly of the free surface height compared to the 14 day mean
ns_maxAmpl = ds_ns.zeta[:].max('ocean_time')-ds_ns.zeta[:].mean('ocean_time')
Max tidal height amplitude¶
plot_max_ampl(ns_maxAmpl,'high res + rough topo',hr_maxAmpl,'high res')
M2 O1 height amplitudes¶
ds_ns = grid_ttide(ds_ns,ds_ns,50)
plot_M2O1(ds_ns,'2x hr + rough topo',ds_hr,'high res')
Scaling tidal height amplitiude of forcing¶
Simple solution would be to scale down the Pamp amplitude of the forcing by simply setting a scaling factor in varinfo.dat. However we have have to think about a justification for that. To determine a rough scaling factor, I compare the ratios between roms and tpxo of sping tide amplitude, M2 and K1.
Determining the scaling factor¶
... so Roms spring tide is 1.55 times to big and needs to be scaled by 1/1.55 = 0.65.
Results:¶
#look at results
file_path = os.path.join(os.environ.get('rawdir'),'waom5_scale','ocean_his_0002.nc')
ds_sc = xr.open_mfdataset(file_path,chunks={'xi_rho':200,'eta_rho':200})
# calculate maximum anomaly of the free surface height compared to the 14 day mean
sc_maxAmpl = ds_sc.zeta[:].max('ocean_time')-ds_sc.zeta[:].mean('ocean_time')
plot_max_ampl(sc_maxAmpl,'high res and scaled forcing',hr_maxAmpl,'base case')
M2 O1 height amplitudes¶
ds_sc = grid_ttide(ds_sc,ds_sc,50)
plot_M2O1(ds_sc,'2x hr + scaling the forcing ampl',ds_hr,'high res')
Mask forcing south of 60S¶
Our forcing is already the reaction of a model to astronomical forcing. The kelvin waves on the shelf might be forced two times in our model. First by the deep ocean wave approaching the shelf and second by the prescribed tides on the shelf. By masking everything south of 60S and ramping full foring to 50S we hope to avoid the two sources of shelf tides.
Generating the mask¶
#mask tidal forcing
import xarray as xr
import os
import sys
import matplotlib.pyplot as plt
pott_path = os.path.join(os.environ.get('rawdir'),'waom10_tidal_mask','Data','waom10','waom10_ptds_large.nc')
pott_ds = xr.open_dataset(pott_path)
grid_path = os.path.join(os.environ.get('rawdir'),'waom10_tidal_mask','Data','waom10','waom10_grd_large.nc')
grid_ds = xr.open_dataset(grid_path)
Pamp_tmp = pott_ds.tide_Pamp.where(grid_ds.lat_rho > -65 ,0.0)
scale = (grid_ds.lat_rho + 65.0)/10
Pamp_masked = Pamp_tmp.where((grid_ds.lat_rho < -65)|(grid_ds.lat_rho > -55),Pamp_tmp*scale)
plt.close()
Pamp_masked[0].plot(vmax=0.5)
plt.show()
Results¶
#look at results
file_path = os.path.join(os.environ.get('rawdir'),'waom10_tidal_mask','ocean_his_tm_65_0006.nc')
ds_tm_65 = xr.open_mfdataset(file_path,chunks={'eta_rho':100,'xi_rho':100})
# calculate maximum anomaly of the free surface height compared to the 14 day mean
tm_65_maxAmpl = ds_tm_65.zeta[336:].max('ocean_time')-ds_tm_65.zeta[336:].mean('ocean_time')
plot_max_ampl(tm_65_maxAmpl,'mask forcing south of 65S',tm_maxAmpl,'mask forcing under ice')
ds_tm_65 = grid_ttide(ds_tm_65,ds_tm_65,50)
plot_M2O1(ds_tm_65,'mask forcing north of 65S',ds_tm,'mask forcing under ice')
... M2 seems to be reuced big time, however O1 is still to big where the coastline is rough. Next try just masking onshelf areas with ramping to deep.
Pamp_tmp = pott_ds.tide_Pamp.where((grid_ds.h > 1000),0.0)
Pamp_tmp = Pamp_tmp.where((grid_ds.zice==0.0),0.0)
scale = (grid_ds.h - 1000)/3000
Pamp_masked = Pamp_tmp.where((grid_ds.h < 1000 )|(grid_ds.h > 4000),Pamp_tmp*scale)
plt.close()
Pamp_masked[4].plot(vmax=0.5)
plt.show()
Results¶
#look at results
file_path = os.path.join(os.environ.get('rawdir'),'waom10_tidal_mask','ocean_his_tm_shelf_0006.nc')
ds_tm_shelf = xr.open_mfdataset(file_path,chunks={'eta_rho':100,'xi_rho':100})
# calculate maximum anomaly of the free surface height compared to the 14 day mean
tm_shelf_maxAmpl = ds_tm_shelf.zeta[336:].max('ocean_time')-ds_tm_shelf.zeta[336:].mean('ocean_time')
plot_max_ampl(tm_shelf_maxAmpl,'mask forcing on shelf',tm_maxAmpl,'mask forcing under ice')
... This looks already a lot better. Any following case will include some sort of shelf mask.
ds_tm_shelf = grid_ttide(ds_tm_shelf,ds_tm_shelf)
plot_M2O1_diff(ds_tm_shelf,'mask on shelf forcing',dsf,'TPXO',vmin=-0.50,vmax=0.50)
Smoothing the mask edge¶
... Just a bit smoother transition at the forcingn/no-forcing edge.
from scipy.ndimage.filters import gaussian_filter, uniform_filter
pott_path = os.path.join(os.environ.get('prodir'),'waom10_tds.nc')
pott_ds = xr.open_dataset(pott_path)
grid_path = os.path.join(os.environ.get('prodir'),'waom10_grd.nc')
grid_ds = xr.open_dataset(grid_path)
mask_tmp = grid_ds.mask_rho.where((grid_ds.zice == 0),0.0)
mask_tmp = mask_tmp.where(grid_ds.h > 1000,0.0)
mask_tmp.values = gaussian_filter(mask_tmp,20)
plt.close()
mask_tmp.plot(size=5)
plt.show()
out_path = os.path.join(pott_path,os.pardir,'waom10_tds_mask.nc')
pott_ds['tide_Pamp']=pott_ds.tide_Pamp*mask_tmp
pott_ds.to_netcdf(out_path,'w')
Results¶
#look at results
file_path = os.path.join(os.environ.get('rawdir'),'waom10_tidal_mask','ocean_his_tm_smooth_0006.nc')
ds_tm_smooth = xr.open_mfdataset(file_path,chunks={'xi_rho':100,'eta_rho':100})
ds_tm_smooth =grid_ttide(ds_tm_smooth,ds_tm_smooth,50)
plot_M2O1_diff(ds_tm_smooth,'mask on shelf smooth',dsf,'TPXO',vmin=-0.50,vmax=0.50)
Drag revisited: tidal mask + 10xCd¶
We apply the same smooth on shelf mask, but increase the bottom drag coefficent to 10 times its normal value.
CAUTION: this means 20x Cd under ice shelves. CAUTION: this is most likely mess up my basal mass balance.
#load roms output that has been forced with this tpxo file
file_path = os.path.join(os.environ.get('rawdir'),'waom10_Cd','ocean_his_x10Cd_0006.nc')
Cd = xr.open_mfdataset(file_path,chunks={'xi_rho':100,'eta_rho':100})
ds_Cd = grid_ttide(Cd,Cd)
plot_M2O1_diff(ds_Cd,'10xCd',dsf,'TPXO',vmin=-0.50,vmax=0.50)
swf bug fix: tidal mask + increased stratification¶
In addition to the tidal mask smooth on shelf, we corrected a bug, which results basically in increased stratification on the shelf.
#look at results
file_path = os.path.join(os.environ.get('rawdir'),'waom10_tidal_signal','ocean_his_hourly.nc')
ds_swf = xr.open_mfdataset(file_path,chunks={'xi_rho':100,'eta_rho':100}).isel(ocean_time=np.arange(365,721))
ds_swf = grid_ttide(ds_swf,ds_swf)
plot_M2O1_diff(ds_swf,'mask on shelf smooth + more stratification',dsf,'TPXO',vmin=-0.50,vmax=0.50)
longer timeseries: tidal mask + increased stratification 30 days¶
Same case as above, just with a 30 days timeseries instead of 14 days as all other cases.
#look at results
file_path = os.path.join(os.environ.get('rawdir'),'waom10_tidal_signal','ocean_his_hourly.nc')
ds_swf_long = xr.open_mfdataset(file_path,chunks={'xi_rho':100,'eta_rho':100})
ds_swf_long = grid_ttide(ds_swf_long,ds_swf_long)
plot_M2O1_diff(ds_swf_long,'as last, but 30 days',dsf,'TPXO',vmin=-0.50,vmax=0.50)
Adjusting the mask edge location¶
To further reduce the jump in O1 amplitude at the mask edge, I push the mask further off shelf.
from scipy.ndimage.filters import gaussian_filter, uniform_filter
pott_path = os.path.join(os.environ.get('rawdir'),'waom10_tidal_mask','Data','waom10','waom10_ptds_large.nc')
pott_ds = xr.open_dataset(pott_path)
grid_path = os.path.join(os.environ.get('rawdir'),'waom10_tidal_mask','Data','waom10','waom10_grd_large.nc')
grid_ds = xr.open_dataset(grid_path)
mask_tmp = grid_ds.mask_rho.where((grid_ds.zice == 0),0.0)
mask_tmp = mask_tmp.where(grid_ds.h > 3000,0.0)
mask_tmp.values = gaussian_filter(mask_tmp,30)
plt.close()
mask_tmp.plot(size=5)
plt.show()
out_path = os.path.join(pott_path,os.pardir,'waom10_ptds_tm_smooth_deep.nc')
pott_ds['tide_Pamp']=pott_ds.tide_Pamp*mask_tmp
pott_ds.to_netcdf(out_path,'w')
Results¶
#look at results
file_path = os.path.join(os.environ.get('rawdir'),'waom10_tidal_mask','ocean_his_tm_smooth_deep_0006.nc')
ds_tm_smooth_deep = xr.open_mfdataset(file_path,chunks={'xi_rho':100,'eta_rho':100})
ds_tm_smooth_deep = grid_ttide(ds_tm_smooth_deep,ds_tm_smooth_deep)
plot_M2O1_diff(ds_tm_smooth_deep,'mask on shelf smooth deep',dsf,'TPXO',vmin=-0.50,vmax=0.50)
Conclusion¶
Nothing seems to permanantly solve the high amplitudes around the shelf area. We go with a mask-on shelf-smooth case, since the other options just increase complexity without much enhancement.